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Abstract 

We investigate a population dynamics model that exhibits a Neimark Sacker bifurcation with a 
period that is naturally close to 4. Beyond the bifurcation, the period becomes soon locked at 4 
CN| ■ due to a strong resonance, and a second attractor of period 2 emerges, which coexists with the 

first attractor over a considerable parameter range. A linear stability analysis and a numerical 
investigation of the second attractor reveal that the bifurcations producing the second attractor 
occur naturally in this type of system. 
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1. Introduction 



The population dynamics of biological species that produce offspring once in their lifetime (semel- 
parous species) provides a generic system that can display a Neimark Sacker bifurcation. A re- 
markable example is the pacific sockcyc salmon (Oncorhynchus nerka), which returns from the 
ocean to spawn in its native lake and dies afterwards. The sockeye fry spend one season in the 
lake, feeding on zooplankton and being eaten by predator fish such as rainbow trout, and then 
migrate to the ocean, from where they return to spawn at the age of 3, 4, or 5 years. Several 

fj^ , sockeye populations display large-amplitude oscillations with a period corresponding to the domi- 

t^J" ' nant generation time of these fish, a phenomenon generally referred to as cyclic dominance [l| . In 

contrast to conventional predator-prey systems, the dynamics of the sockeyc-trout system in the 

f**^ ■ lake is piecewise continuous in time, due to the salmon migrating to the ocean and returning at 

the age of 4 (±1) years. This mechanism, which puts the model in the context of delay dynamical 
systems, will turn out to be one of the crucial ingredients. The second one is the external periodic 
drive caused by the yearly rhythm of seasons. The ratio of these two time scales is thus fixed by 
an internal mechanism which turns out to be vital for the observed phenomena. 

r% \ While conventional predator-prey systems typically display a Hopf bifurcation |2j , our piecewise 

continuous system displays a Neimark Sacker bifurcation [3j, i.e., an oscillatory instability in ex- 
ternally driven systems. Since the dominant life span of the salmon is 4 years, the period of the 
oscillation is naturally close to 4, and when the control parameter is further increased beyond the 
bifurcation, the period becomes locked exactly at 4 due to the strong resonance familiar from the 
theory of the Neimark Sacker bifurcation [4| . 

Resonances and the associated synchronisation phenomena are common features in driven systems, 
with widespread applications in science and technology (cf. [5[ for a recent review). But unlike 
strong resonances one normally requires a fine tuning of parameters, the subtle synchronisation 
phenomena like phase locking or universal scaling behaviour commonly related with the quasipcri- 
odic transition to chaos [6|,|7| may be sensitive to perturbations. In contrast, strong resonances in 
Neimark Sacker bifurcations generate very stable periodic motion over large regions in parameter 
space, and thus can provide a vital mechanism for cyclic dominance in biological systems where 
noise and parameter fluctuations are prevalent. 

In this paper, we will explore further a model for the population dynamics of the pacific sockeye 
salmon. We will show that it can display a second attractor, which has the period 2 and coexists 
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with the attractor of period 4. Due to the strong resonance, the second attractor is clearly visible 
in bifurcation diagrams where only one of the 6 dimensions of the system is plotted as function of 
a control parameter. By using computer simulations, we will show that the second attractor arises 
through a flip bifurcation followed by a subcritical Ncimark Sackcr bifurcation, and by using a 
theory that is linearised around the fixed point as well as around the attractor of period 4, we will 
show that the flip bifurcation can naturally occur soon after the first Neimark Sacker bifurcation. 

2. The model 

The dynamics of the biomass of sockeye fry s n (t), of their predator p n (t), and of their zooplankton 
food z n (t) in year number n during the growth season from spring (t = 0) to fall (t = T) are given 
by the following set of coupled differential equations [3( : 

— S n (t) = Xg sz (s n (t),Z n (t)) - g P s(Pn(t),S n (t)) ~ d a ■ S n (t) 

j t z n (t) = z n (t)-(l-^-j -g sz (s n (t),z n (t)) (1) 

— p n (t) = Xg pa (jp n (t), S n (t)) - dp ■ p n (t) 

The terms in these equations represent the biomass changes due to eating and being eaten, and 
the losses due to respiration and death. The two feeding terms are given by 



9sz(s,z) = 
9 P s(p,s) = 


1 + c s ■ s + z 

s ■ p 


ps 1 + c p ■ p + s 



(2) 

These functional forms include saturation at high prey densities, and a predator interference term 
in the denominator (Bcddington functional response [slSl]). These equations of motion determining 
the time evolution of the biomass during a single season are supplemented with matching conditions 
used to determine the biomasses of the species at the beginning of the next season from their values 
at the end of the previous season(s). To keep the model as simple as possible we assume that no 
major change in the population of predators occurs. The zooplankton level at the end of one year 
has no effect on the following year lfjj, but its saturation value in the next season, K n +i, shows a 
dependence on the nutrient input due to the number of spawning adults, and thereby on s n +i(0) 
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z n+1 (0) = K n+1 =K +(s S : +li0 \ n , ) (3) 

Pn+l(0) = Pn (T) (4) 

Most importantly the matching condition for the sockeye fry population contains the proportion 
ei and e 2 of surviving sockeye that return to their native lakes at the age of 3 and 5 to spawn and 
die: 

s«+i(0) = 7[(1 - ei - e 2 )s„_3(T') + eiS n _ 2 (T) + e 2 s n - 4 (T)} (5) 

On the one hand the model may be considered as a very simple example of a pieccwisc smooth 



dynamical system |13| . where the discontinuities are contained in the matching condition between 
seasons. On the other hand the model may be considered as a periodically driven non autonomous 
system where the external drive is provided by the seasons. It is therefore an obvious approach to 
consider the stroboscopic map. If we denote by 3>[o,t] the fi° w of the differential equation (fTJ) the 
biomasses at the beginning and at the end of the year n + 1 are related by 

(s n+ i(T),p n+1 (T),z n+ i(T)) = $ [o ,T](sn+i(0),p„+i(0),^ n+ i(0)) (6) 

The condition ([3]) tells us that the zooplankton is entirely determined by sockeye population so that 
z n drops from Eq.(|6]) as a dynamical variable. Taking Eqs.Q and ([5]) into account the remaining 
biomass evolution can be written as 

s„ + i(T) = h s {e 1 s n - 2 {T) + (1 - ei - e 2 )s n ^ 3 (T) + e 2 s n -i{T),p n (T)) 

p n +i(T) = h p (eis n -2(T) + (1 - ei - e 2 )s n - 3 (T) + e 2 s„-4,(T),p n (T)) . (7) 



Thus we arrive at a six dimensional time discrete dynamical system. The computation of the right 
hand sides, h s and h p , requires the integration of the equations of motion, Eq.([T]), a task which can 
be accomplished easily by numerical means. However, the structure of the map ([7]), induced by the 
matching condition ([5]) will turn out to be crucial for the properties of the population dynamics. 

3. The main bifurcation generating cyclic dominance 

It is our main purpose to pinpoint an underlying generic mechanism which creates stable cyclic 
motion among the populations in our model. For that purpose we employ a brief bifurcation 
analysis for realistic combinations of the parameters of the model. Given the fairly large number 
of parameters in our model and the fairly large number of degrees of freedom we just resort to 
plain numerical simulations. In fact, employing for our specific purpose numerical continuation 
tools like matcont |14| our auto [15[ would be using a sledge-hammer to crack a nut. In addition, 
our specific numerical study mimics approaches which could be used in evaluating empirical data, 
and thus could be a guide for their analysis. 

Figure Q] shows the bifurcation diagram for a realistic set of parameter values. For each value 
of the parameter <5, several initial conditions were drawn at random from uniform probability 
distributions on a 6-dimcnsional cube (0; 10]. When coexisting attractors were found, these were 
numerically continued until they became unstable. After a transient time of 5000 years (longer close 
to bifurcations), the values for s n (T) were plotted for 20 additional years. As 6 increases, the fixed 
point becomes unstable due to a Neimark Sacker bifurcation, and a quasiperiodic attractor arises. 
At S ~ 4.47, trajectories can become locked exactly at the period 4, but the quasiperiodic attractor 
remains stable until 5 ~ 4.6. This picture is consistent with a Neimark Sacker bifurcation close to 
a strong resonance of order 1:4, which is known to display quite a subtle bifurcation structure [4|, 
including the creation of homoclinic tangles. 




Figure 1: Bifurcation diagram of the biomass of the sockeye fry at the end of their first year. Parameter values used 
for this diagram are: A = 0.85, a sz = 10, a ps = 0.2, c s = 1, e p = 0.2, d s = 1.1, dp = 0.3, 7 = 0.2, e\ = 0, £2 = 0.1, 
K = 2.5, and 5 = 2. 



At 8 ~ 5.36, a second attractor with period 2 appears. Such a feature is normally not related 
with the generic strong resonant Neimark Sacker bifurcation mentioned above. The two arms 
of this attractor appear as if they are created as an unstable period two orbit at S ~ 5.1 via 
a flip bifurcation from the already unstable fixed point. The unstable arms then become stable 
at 6 ~ 5.36 when the complex conjugate pair of eigenvalues moved back into the unit circle. 
Even without continuation tools tracking saddle orbits we are able to vindicate this picture by 
evaluating the size of the basin of attraction of the period two attractor and by monitoring the 
transient dynamics. 

Figure [2] shows an example of the transient dynamics in the stable period two regime. Since only 
every other iteration is displayed the transient consists of a damped oscillation with a period close 
to 4. By measuring the decay rate r of the oscillations around the stable period two orbit for 
different values of S we determined the bifurcation value 5 C with an accuracy of 10~ 6 . Sufficiently 
close to the bifurcation, r goes linearly to as the instability is approached (lower inset in Figure 



[2|) . We then estimated the diameter of the basin of attraction by using a one dimensional cross 
section along the direction of the first salmon population. Initial conditions were distributed along 
this cross section, and those trajectories which return to the period two state determine the extent 
of the basin along this direction in phase space. Figure [3] shows the result for the size of the 
basin, and the characteristic square root law with regard to the distance from the bifurcation point 
is recorded. Thus the bifurcation which causes the stability of the period two state is clearly a 
subcritical bifurcation. 
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Figure 2: Trajectory approaching the period 2 attractor. <5 = 5.45, all other parameters are as in figure 1. Only 
every second data point is plotted. Upper inset: close-up of the trajectory. Lower inset: decay rate of the oscillations 
when approaching the instability. The dashed line is a linear function fitted to the numerically determined decay 
rates, yielding 5 C = 5.356111 ± 5 ■ 10" 7 . 

All these results demonstrate that the second attractor emerges due to a subcritical Neimark Sacker 
bifurcation, and this secondary Neimark Sacker bifurcation is intimately related to the bifurcation 
creating the stable period four attractor. We are now going to uncover the underlying mechanism 
by analytical considerations. 
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Figure 3: Size of the basin of attraction of the period 2 attractor. Parameters are as in figure[T] and <5 C = 5.356111. 
The dashed line is a power law with exponent 1/2. 



4. Linear stability analysis 

A linear stability analysis of the dynamical system close to the first Neimark Sacker bifurcation 
reveals that the scenario outlined in the previous section is indeed a natural scenario for this type 



of system. The continuous population dynamics during the season, together with the matching 
conditions applied between two seasons, can be viewed as a Poincare map, Eq.([7]), giving the 
sockcye and predator biomasscs at the end of one year as function of the biomasses at the end of 
the previous years. Close to the fixed point (s*,p*), the dynamics can be approximated by linear 
terms. Denoting the distance of the biomasses from their fixed point value by Ss n = s n (T) — s* 
and Sp n = p n (T) — p*, the linear approximation of Eq.© has the form 



(8) 
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(9) 



denote the partial derivatives of the Poincar map at the fixed point. These parameter values can 
in principle be obtained by a numerical evaluation of the population dynamics, Eq.{T]), close to 
the fixed point. Any consistent underlying dynamical model should give a negative value of the 
parameter /3, since an increase in trout biomass leads to a decrease in salmon biomass. The other 
parameters in Eq.([9]) have to be non- negative, since more parents imply more offspring and since 
an increase in salmon biomass leads to an increase in trout biomass. However, it is notoriously 
difficult to derive such conditions from equations of motion like Eq. ((1|) . 

The eigenvalues of the matrix in Eq. (J8j) determine the nature of the dynamics near the bifurcation. 
When all eigenvalues have an absolute value smaller than 1, the fixed point is stable, and the 
dynamics converges to this fixed point. When the absolute value of one or more eigenvalues is 
larger than 1, the fixed point is unstable, and the dynamics approaches a different attractor. The 
eigenvalues \x satisfy the characteristic equation 



p 5 (P -p) + [eut 2 + (1 - ei - e 2 )p + e 2 ][/3a - a{fi - n)] = . 



(10) 



In order to understand the possible bifurcation scenarios, we first consider the case that only the 
parameters a, (3 and j3 are nonzero. This means that the four salmon lines are independent from 
each other, ei = e 2 = 0, and from the trout, a = 0. The eigenvalues of the matrix are j3, 0, and 
the four fourth roots of a. The eigenvalue j3 characterises the time evolution of the trout, which is 
decoupled from the salmon in the considered case. Global stability criteria require this eigenvalue 
to be smaller than one, < {) < 1. The eigenvalue is due to the fifth column of the matrix being 
redundant if e 2 = 0. The four fourth roots of a describe the salmon dynamics in the vicinity of 
the fixed point. The sequence Ss n has trivially the period 4 and simply iterates the initial four 
values, with an amplitude decreasing for a < 1 and increasing otherwise. When a is increased 
from a value smaller than 1 to a value larger than 1, all four eigenvalues a 1 ' 4 cross the unit circle 
simultaneously, and the fixed point becomes unstable. This degeneracy is lifted when the other 
parameters, e±, e 2 , and a, are made nonzero. As long as these parameters are not large, one can 
expect the four main eigenvalues to remain close to the unit circle, implying a (possibly damped) 
oscillation with a period close to 4. If the complex conjugate pair of eigenvalues is the first one 
to cross the unit circle, a Neimark Sacker bifurcation occurs. If the negative eigenvalues crosses 
the unit circle first, a flip bifurcation occurs. For the parameter values chosen in the simulation 
(Figure [1]), the Neimark Sacker bifurcation is observed first. 

The entire bifurcation scenario can be computed from Eq. (|10j) using the Schur-Cohn-Jury criterion 
[16j , but the main mechanism which is at the heart of the bifurcation discussed in the previous 
section can be already uncovered if we resort to a perturbative approach. In order to gain some 
analytical insights, we assume that the parameters ei, e 2 , and a are so small that the change in 
\x can be calculated approximately by performing a first order Taylor expansion in a — 1, close 
to the bifurcation. By implicit differentiation of Eq. (fTOjl we obtain for the relative change of the 



eigenvalues close to the unit circle 
dfi a — I 



/'o 



4(1-/3) 
ei +e 2 


2 
ei +e 2 



(-/3)a 



4(1 + /?! 



if /j = 1 

if Mo = -1 (11) 



(-/3)a/3 i / (-/3)q 

T t I ei ~ £2 + ,,_, ; 59 , if Mo = ±* 



4(1 + /3 2 ) 4 V 4(1 + /3 2 



where //o denotes the eigenvalue of the unperturbed case. The real part of this ratio d/i/no de- 
termines the stability properties of the eigenmode. As long as the parameter /3 is not too small, 
i.e. as long as the decay of the trout population is not massive, the negative contribution to the 
positive real eigenvalue (the case /Iq = 1) dominates, and this mode remains stable on changing 
the bifurcation parameter a — 1. Thus no saddle node bifurcation is expected to occur. The order 
of the other two cases, a bifurcation caused by a negative eigenvalue (the case /l*o = ~ 1) an d a 
strongly resonant Ncimark Sacker bifurcation (the case /j,q = ±i) is now determined by the balance 
between the unfolding parameters. If ei +e 2 is sufficiently large compared to a, the Neimark Sacker 
bifurcation appears first, as the negative contribution to the second case in Eq. (fTTj) is larger (in 
modulus) than the negative contribution to the third case. In quantitative terms Eq. (|ll[) yields 
for the suppression of the saddle node bifurcation the inequality 

Y^ >£l + e2 (12) 

whereas the Neimark Sacker bifurcation precedes the flip bifurcation if 

R3)a(l-i) 
61 +£2> (1 + #U + ^) (13) 

is satisfied (remember that /3 < and < /3 < 1). 

These considerations show that as long as the coupling between the different salmon lines and 
between the salmon and the trout is small compared to the coupling between salmon parents and 
offspring (i.e. e\ <C 1, e 2 <C 1, and a <C 1), the oscillation period of the population is close to four. 
If the coupling between the salmon lines is finite and beyond a threshold value determined by the 
small coupling between the salmon and the trout (i.e. by the parameter a, cf. Eq. (|13[) ) the period 
four oscillating mode becomes unstable first as the bifurcation parameter a is increased, and the 
complex conjugate eigenvalues close to the imaginary axis arc pushed towards the unit circle. The 
complex conjugate eigenvalues are therefore the first ones to cross the unit circle followed by the 
eigenvalue on the negative real axis. But we require as well a decay rate for trouts which is not 
excessively large, Eq.(fT2"|). to prevent the saddle node bifurcation to occur. In fact, in the limit 
(3 — >• the two conditions, Eqs.(fT2"|) and (fl~3|) become mutually exclusive. Thus the dynamics of the 
predators play an important role and the scenario would not appear if an adiabatic elimination of 
the trout population would have been an option. The scenario found in the computer simulations, 
where first a Neimark Sacker bifurcation occurs and then a flip bifurcation, is thus a generic scenario 
in this type of system. Clearly, our linear analysis cannot keep track of the stability properties 
of the period two orbit created in the secondary flip bifurcation. Certainly this orbit inherits 
its stability properties from the bifurcating fixed point, but we cannot predict without a tedious 
normal form calculation whether the two unstable complex conjugated eigenvalues move back into 
the unit circle, although such a re-entrance phenomenon is not uncommon. While we observed the 
subcritical Ncimark Sacker bifurcation on the period-2 orbit for a broad set of parameter values, it 
did not occur in other regions of parameter space. But the presented analytical considerations have 
shown that the stable period two and period four oscillations arc caused by the same degenerate 
Neimark Sacker bifurcation, and we were able as well to derive criteria for the occurrence of this 
scenario. 

5. Conclusions 

We have shown that a three-species population dynamics model can serve as a model system for 
observing the coexistence of two attractors, one of period four and one of period two. This is 
remarkable, since this type of scenario is of codimension 3, requiring a Neimark-Sacker bifurcation 



occurring with a period close to four, and a third eigenvalue crossing the unit cycle for a nearby 
value of the control parameter. We have identified conditions such that the bifurcation of high 
codimension does not require fine tuning of additional parameters. First of all, the external stimulus 
of yearly season and the return of species for spawning provides the strong resonance condition of 
the Neimark Sacker bifurcation and the reduction to a time discrete dynamical model. The second 
important ingredient is a predator population which actively takes part in the dynamics and 
cannot be taken into account on a plain adiabatic description. Finally a weak coupling between 
the different spawning populations is required. While the zooplankton provides the source for 
keeping the nonequilibrium dynamics going, it does not independently take part in the dynamics 
but depends solely on the current salmon population. While our considerations have been based 
on a simple but realistic model, Eq.([T]), we should emphasise that the conclusions have a fairly 
broad remit as the time discrete dynamics (0 only takes a few generic features of the model into 
account. Population models that have discrete generations can thus represent a natural situation 
for a scenario that was previously expected to require very improbable special conditions. 
At a more fundamental abstract level the model proposed in this contribution is at the interface 
between piecewise smooth dynamical systems and driven delay dynamics, both being fields of 



active research in applied dynamical systems theory [13j, lUj . Although being a standard aspect 
of bifurcation theory, the mechanism presented here could have a wider impact for the study of 
robust periodic motion and synchronisation in a large class of experimentally relevant setups in 
science and engineering. 
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